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Abstract 

Polarized protons have never been accelerated to more than about 25GeV. 
To achieve polarized proton beams in RHIC (250GeV), HERA (820GeV), and 
the TEVATRON (900GeV), ideas and techniques new to accelerator physics are 
needed. In this publication we will stress an important aspect of very high en- 
ergy polarized proton beams, namely the fact that the equilibrium polarization 
direction can vary substantially across the beam in the interaction region of a 
high energy experiment when no countermeasure is taken. Such a divergence 
of the polarization direction would not only diminish the average polarization 
available to the particle physics experiment, but it would also make the polar- 
ization involved in each collision analyzed in a detector strongly dependent on 
the phase space position of the interacting particle. In order to analyze and 
compensate this effect, methods for computing the equilibrium polarization di- 
rection are needed. In this paper we introduce the method of stroboscopic 
averaging, which computes this direction in a very efficient way. Since only 
tracking data is needed, our method can be implemented easily in existing spin 
tracking programs. Several examples demonstrate the importance of the spin 
divergence and the applicability of stroboscopic averaging. 
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Introduction 

In order to maximize the number of collisions of stored particles in a storage ring 
system one tries to maximize the total number of particles in the bunches and tries to 
minimize the emittances so that the particle distribution across phase space is narrow 
and the phase space density is high. At equilibrium the phase space distribution does 
not change in time and is therefore periodic in the machine azimuth. 

If, in addition, the beam is spin polarized, one requires that the polarization is 
high. As first emphasized by D. Barber for energies of the order of ITeV, a 
fundamental limitation to the polarization of particle beams becomes important. To 
put our work in context we repeat the arguments here. 

Spins traveling with particles in electromagnetic fields precess according to the 
Thomas-Bargmann-Michel-Telegdi equation (T-BMT) to be discussed below. The 
guide fields in storage rings are produced by dipole and quadrupole magnets. The 
dipole fields constrain the particles to almost circular orbits and the quadrupole fields 
focus the beam, thus ensuring that the particles do not drift too far away from the 
central orbit. 

In horizontal dipoles, spins precess only around the vertical field direction. The 
quadrupoles have vertical and horizontal fields and additionally cause the spins to 
precess away from the vertical direction. The strength of the spin precession and the 
precession axis in machine magnets depends on the trajectory and the energy of the 
particle. Thus in one turn around the ring the effective precession axis can deviate 
from the vertical and will depend on the initial position of the particle in six dimen- 
sional phase space. From this it is clear that if an equilibrium spin distribution exists, 
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i.e. if the polarization vector at every phase space point is periodic in the machine 
azimuth, it will vary across the orbital phase space. This field of equilibrium spin 
directions in phase space does not change from turn to turn when particles propagate 
through the accelerator, although after each turn the particles find themselves at 
new positions in phase space. These directions, which we denote by the unit vector 
n{z, 6), where z denotes the position in the six dimensional phase space of the beam 
and 9 is the generalized azimuth, was first introduced by Derbenev and Kondratenko 
0] in the theory of radiative electron polarization. Note that n{z, 6) is usually not 
an eigenvector of the spin transfer matrix at some phase space point since the spin 
of a particle changes after one turn around the ring, but the eigenvector would not 
change. 

Thus once we know this direction n{z,6), the phase space dependent polariza- 
tion p{z, 9) in this direction, and the phase space density function p{z^ 9) we have a 
complete specification of the polarization state of a beam of spin 1/2 particles. 

The maximizing of the polarization of the ensemble implies two conditions; the 
polarization p{z, 9) at each point in phase space should be high and the polarization 
vector n(z, 9) at each point should be almost parallel to the average polarization 
vector of the beam. 

According to the T-BMT equation, the rate of spin precession is roughly propor- 
tional to 07 where a = ((7 — 2)/2 is the anomalous part of the spin g factor and 7 is 
the Lorentz factor. At very high energy, as for example in the HERA proton ring [0, 
it could happen that on average n{z, 9) deviates by tens of degrees from the phase 
space average of n. Thus even if each point in phase space were 100% polarized the 
average polarization could be much smaller than 100%. Clearly it is very important 
to have accurate and efficient methods for calculating n(z, 9) and for ensuring that 
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the spread of n(i', 6) is as small as possible. 

However, although it has been straightforward to define n{z, 9) , this vector is 
not easy to calculate in general and much effort has been expended on this topic, 
but mainly for electrons at energies up to 46GeV. Except for the Fourier expansion 
formalism introduced in all other methods developed so far are explicitly pertur- 
bative, and either do not go to high enough order 0, |^, |^ or have problems with 
convergence at high order and high energy |0, H. 

In this paper we describe a new method for obtaining n{z,6). It is based on 
multi-turn tracking and the averaging of the spin viewed stroboscopically from turn 
to turn. Since this innovative approach only requires tracking data, it is fast and 
very easy to implement in existing tracking codes. We will show that the convergence 
speed promises rather quick execution when simulating realistic accelerators. How- 
ever, probably the main advantage over other methods is the fact that stroboscopic 
averaging does not have an inherent problem with either orbit or spin orbit reso- 
nances due to its non-perturbative nature. This allows the behavior of the periodic 
spin solution close to resonances to be analyzed. 

1 The Spin— Orbit System 

The motion of the spin of particles traveling in electromagnetic fields is governed by 
the equations of motion 



-{6) = Q{m,d)xm- (2) 



Here 6 is an independent variable parametrizing the d dimensional particle phase 
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space trajectory z{6) and the spin trajectory s{6). In circular accelerators 6 is the 
azimuth. The rest frame spin vector s has 3 components and we normally deal with 
orbital phase space vectors z which have 6 components. In accelerator physics these 
components are usually the positions and momenta of a particle combined with its 
energy and the time of flight. 

We are neglecting the Stern-Gerlach forces since they are very small in comparison 
with the Lorentz force. Equation (§) is the Thomas-BMT equation along an orbit 
parametrized by z{6) [P, Because we deal with a circular accelerator at fixed 
energy, v{z,9) and Q{z,9), which depend on the guide fields, are periodic in 9 with 
period 27r corresponding to the circumference of the ring. Due to the precession 
equation (0) the length of the spin vector s does not change along the azimuth. 

The dynamical system (|1|J^) allows to formulate the following partial differential 
equation for the evolution of a field f{z, 9): 

— » 

where the three components of the solution / depend on and 9. In our applications 
/ will describe the propagation of a spin distribution associated with a particle beam 
and this physical interpretation can be adopted because of the following. A solution 
/(i*, 9) to equation can be found by specifying an arbitrary f{z, 9q) at initial 
azimuth 6*0 and propagating it to 9 by integrating equations (|^) and (^. In fact 

s{9) = Rm,0) (4) 
solves equation if / solves equation (Q). We say that / is normalized if 



I/I = \//f + /! + /! = 1- (5) 



We call every normalized solution of equation (j^) a 'spin field'. 
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The 'n-axis' introduced in the Introduction is a special spin field which is periodic 
in 9 with period 27r 0: 

n(f, e + 2Ti)= n{z, 9) . (6) 

Since equation (|^) represents a pure rotation the propagation of the spin vector can 
be described by a 3 x 3 orthogonal matrix. We denote this rotation matrix which 
propagates initial spins s{9q) along a given orbit trajectory z{9) by R{z{9q), 9, 9q) so 
that 

m = R{m).eM-s{o^) ■ (7) 

Because for a spin field / the spin trajectory (|) solves equation (0), we get 

Am, 0) = mm), Oo) ■ moo), oo) . (s) 

If f{z, 9) is a spin field, then f{z, 9 + 27rm) is a spin field (m= integer). This follows 
from equation @ because v and fl are periodic in 9. Thus equation (||) generalizes 
to 

f{z{9), + 2Txm) = R{z{9o), ^, ^^o) ■ f{m),0o + 27rm) . (m = integer) (9) 
Since an n-axis is a periodic spin field, we observe by and that 

+ 27r), ^) = E(^(^), 9 + 2n, 9) ■ n{z{9), ^) • (10) 



Alternatively this equation can be used for defining the n-axis ||rT|. The matrix 
R{z{9), 9 + 2tt, 9) is called the one turn spin transfer matrix for the trajectory z{9). 

In the special case where the orbital motion is determined by a Hamiltonian we 
have 

v{z,d) = {z,Horb{z,9)} , (11) 

where Horb denotes the orbital Hamiltonian. Furthermore in this case one can define 
for the whole spin-orbit system a Hamiltonian given by [|l^ 



H{z, s, 9) = Horbiz, 9) + Hsp^n{z, s, 9) , (12) 
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where 

Hspin{z, s,6) = s- Q{z, 6) . (13) 

The Poisson brackets of this Hamiltonian lead to equations (|ip and (0) if the Stern- 
Gerlach forces are nea; lected III. 



2 Construction of Periodic Spin Fields by Strobo- 
scopic Averaging 

To find solutions of equation @ which are periodic in 6 by our new method, one first 
constructs an arbitrary spin field /. One then constructs the following 'stroboscopic 
average' of /: 

< / > (i-, 9) = „lim [-1- f: m + 27rm)] . (14) 

Since the convergence and differentiability of the sequence in (|14D can in general not 
be guaranteed, the limit is only taken formally. The problem of the convergence will 
be analyzed in more detail in section From (0) it follows that < / > is formally 
periodic in 6. Moreover, because (^) is a linear equation and because v and Q are 
periodic in 6 we observe for any spin field / that < / > is formally also a solution of 
equation (|^). Hence we come to a first conclusion: 

• If /* is a spin field, then < /* > is a solution of equation (^) which is periodic in 

e. 

If for / the stroboscopic average < / > vanishes nowhere in the d + 1 dimensional 
space, then we define 

< f>norm = < f>/\<f>\. (15) 

In general < / > is not normalized, but the modulus of < / > is conserved and 
< / >norm IS a spiu field which is periodic in 9. Hence we come to the following 
second conclusion of this section: 
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• If / is a spin field with the property that < / > vanishes nowhere in the d + 1 
dimensional space, then < / >norm has all the properties of an n-axis. 

This result shows that an n-axis can be obtained from a spin field / for which < / > 
vanishes nowhere. In the next section we will derive a tracking algorithm based on 
this. 

One practical choice of / is characterized for all z by 

f{z,9o)=no{0o) , (16) 
where no{9) denotes the so called closed orbit spin axis defined by 

no{e)=R{zc.o.{0),0 + 27r,e)-no{e) , (17) 
where ZcoX^) is the closed orbit. 

3 The Tracking Algorithm for the n-Axis using 
Stroboscopic Averaging 

In this section we develop a tracking algorithm which provides an efficient way to 
evaluate an n-axis at z = zo,6 = Oq. Choosing a spin field / and replacing Oq by 
00 — 27rm in equation (^ we get for every integer m 

fiz{9), 9 + 27rm) = E(^(^o - 27rm), 9, 9q - 2nm) ■ f{z{9o - 27rm), ^o) • (18) 

If we choose an orbit with z{9q) = zq, then inserting this into equation (0) results in 

1 

< / > {zo, 9o) = }^[-—r E ^(^(^0 - 27rm), 9o, 9o - 2nm) ■ f{z{9o - 27rm), ^o)] ■ 

(19) 

Normalization of < / > yields an n-axis at [zq, 9q). 
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To apply the tracking algorithm, the infinite sum involved in the stroboscopic 
average (^) is replaced by a finite sum of + 1 terms so that we approximate 



1 ^ 



< / > (ib, ^o) ~ < / >7V (^0, ^o) = jT—r J2 /(^o, Oo + 27rm) 

+ m=0 



1 ^ 

E ^(^(^0 - 27rm), ^0, - ^nm) ■ /(z(^o - 27rm), 



m=0 

which yields the following approximation of the n-axis 



^(^0, ^o) ~ ^ ^ ,1 ■ (21) 

I < / >^ (2:0,^0)1 



The stroboscopic average < / >Ar in equation (|2^) has a very simple physical interpre- 
tation which illustrates its practical importance. If a particle beam is approximated 
by a phase space density, disregarding its discrete structure, then we can associate a 
spin field f{z, 6q) with the particle beam at the azimuth 6q. If one installs a point like 
'gedanken' polarimeter at a phase space point Zq = z{6q) and azimuth 6q, then this 
polarimeter initially measures /(io,^o)- When the particle beam passes the azimuth 
^0 after one turn around the ring, the polarimeter measures 

R{z{eo - 27r), 00, 00 - 27r) ■ f{z{eo - 2n), 60) = f{zo, Oo + 2it) . 

After the beam has traveled around the storage ring N times and the polarization has 
been measured whenever the beam passed the 'gedanken' polarimeter, one only has 
to average over the different measurements in order to obtain < f >]\j. If the particles 
of a beam are polarized parallel to n{z, 60) at every phase space point, then the spin 
field of the beam is invariant from turn to turn due to the periodicity property in 
equation (P). But in addition, even for beams which are not polarized parallel to n, 
we see that the polarization observed at a phase space point z and azimuth ^0 is still 
parallel to n{z,6o), if one averages over many measurements taken when the beam 
has passed the azimuth 6q. 
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For the special choice f{z, 9q) = no(^o) we can simphfy < / >Ar to 

1 ^ 

<f>N {zo, 9o) = —— J2 - 27rm), ^o, - 27rm) ■ no(^o) • (22) 

+ -'- m=0 

Equations (^) and (|21| ) define an algorithm for obtaining an n-axis. We see that 
the only information needed from tracking is the set of the + 1 phase space points 
^•(^0), z{9q-2'k), z{eo-27iN) and the N matrices E(^(^o-27r), Oq, 60-271), R{z{eo- 
47r), ^^0, ^0 — 47r), . . . ,R{z{6q — 2tiN), 9q, 9q — 2tiN). Each matrix is a product of one 
turn spin transfer matrices R{z,9q + 27i,9q). This means that one tracks along the 
orbit z{9) to obtain the spin transfer matrix R{z{9),9o,9) and stores it at the 
instants, where 9 = 9q ~ 2t^N, . . . ,9 = 9q — Att, 9 = 9q — 27r. The function f{z, 9q) 
is chosen independently of the tracking results (for example one can take the choice 
f{z,9o) = no(^o) of equation (|16D and (|2|)). 

The following two kinds of pathologies can occur: 

• The n-axis is not unique: if the proposed algorithm converges, then the result 
could depend on the choice of f{z, 9o). 

• The stroboscopic average < / > vanishes for +00 or the sequence in 
equation (|^) does not converge. 

Both pathologies can be studied with the algorithm. The first situation occurs for 



systems on spin-orbit resonances |Tj]. In all examples studied so far, the strobo- 
scopic average seems to converge, implying the existence of an n-axis. In the second 
situation, the point like polarimeter at zq mentioned above monitors an averaged 
polarization which either vanishes or fiuctuates indefinitely. 



11 



4 Efficient Implementation Only Using One Turn 
Information 



In the previously outlined formalism for evaluating an n-axis by stroboscopic aver- 
aging it became apparent that only knowledge about one turn spin transfer matrices 
is required. One can therefore reformulate the algorithm of section |^ in terms of one 
turn maps which are to be taken at a fixed but arbitrary azimuth value and thereby 
obtain a more practical algorithm. Thus we introduce the one turn orbit transfer map 
M which maps initial coordinates Zi into final coordinates zj = M{zi). Then in our 
notation we have z{Oq + 2tx) = M{z{9q)). To describe the transport of particles with 
spins s, we write for simplification the one turn spin transfer matrix R{z, 9q + 27r, ^^o) 
as R{z) so that we have Sf = R{zi) ■ Si. All other quantities which depend on 9 are 
taken at the specified azimuth 9q. For simplification, this azimuth is not indicated in 
the following. As already mentioned, equation ( jlOD can be used to define n-axes at 
6*0. This condition will be called the periodicity condition and now reads 



4.1 Recipe 

To illustrate the process of evaluating an n-axis at zq and 9q in the case of linear orbit 
motion, we establish a recipe. 

1. Compute the linearized one turn phase space transfer map Zf = M ■ Zi. 

2. Define the set of + 1 phase space points 




(23) 




(24) 
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3. Compute the rotation matrix R{zc,o.) which describes the spin motion for par- 
ticles on the closed orbit Zc.oX^) extract the corresponding rotation vector 
no- This is the periodic spin solution for particles on the closed orbit. 

4. Starting with a spin parallel to at every phase space point in C, track until 
the phase space point zq is reached. For a given j this requires tracking j turns 
starting at cj. 

5. Define the set of spin tracking results as 

B = {bo{zo) = no, b.izo) = R{ci) ■ . . . ■ R{c,) ■ no|j e {1, . . . , iV}} . (25) 

6. Define the sum of the elements in B as sn{zo) = YljLo ^i(^o) and for \sn\ 7^ 
define Hn = sn/\sn\- 

The vector sn{zo) is equivalent to < / >Ar {zq, 6q) in section if the initial distribu- 
tion of spins is given by uq as in equation (|2^). 

4.2 Convergence Speed 

It will now be shown that if the angle between r?o and hj{zo) is smaller than some 
positive number ^ < 7r/2 for all j G {!,..., + 1}, then satisfies the peri- 
odicity condition (|23| ) for the ra-axis up to an error which is smaller or equal to 
2 sec(,^/2) tan(^)/(A^ + 1). Since evaluating B by the recipe of section requires 
tracking T = [N + l)N/2 turns, the accuracy is bounded by y^2/Tsec(,^/2) tan(,^). 
This slow convergence with the square root of T is a very serious limitation and in the 
next section we will demonstrate how the convergence can be considerably improved. 

The proof of this convergence property goes along the following lines. The average 
s/v has been defined by 

I N j 

M^o) = E n Ri^k) ■ no . (26) 
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Here we adopt the convention nfc=i^(cfc) = nfe=i^(cfc_i) = 1 and Y[i=iR{ck) is 
taken to mean the following order of multiplication: R{ci) ■ . . . ■ R{cj). To check how 
well sV satisfies the periodicity conditions (^) of the n-axis we calculate 

1 A' j 



ij ■ no 

j=0 k=l 



I N-l j 

T7-t(^o+ E ER{ck)-no) , (27) 



I N j 

R{zo)-sn{So) = T7-TEn^(cfc)-^o , (28) 

R{zo) ■ sn{zq) - sn{K ■ zq) = ^ ^ ^ (Rizo) ■ bN{zo) - no) 

= ]^(W(M-^o)-^o)- (29) 




Figure 1: Estimation of convergence speed 

The length |6/v+i(il£' -^o) ~no\ is smaller than 2sin(^/2) as shown in figure |l[ The 

length of sjv is at least cos(^); and here it becomes essential that there is limit of 7r/2 
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on the angle ^. This information is sufficient to estabhsh the following inequality 

\R{Zo) ■ UNiZo) - UNiK- Zo)\ = \R{Zo) ■ - ^ 

\sn{zo)\ \sn{M- Zo)\ 
■\R{zo) ■ sn{zo) - sn{M_ ■ Zq) 



\sn\zq) 





\sn{M-zq)\ - 


\sn{zo)\ 


1 


\sn{M-zo)\ 



sn{K- Zo)\ 
\sn{zo)\ 

+ \\R{^o)sn{zo)\ - \sn{M-zo)\\) 
2 

\sn[zo)\ 
< 4sin(e/2) 

- (N + l) cos(0 ■ ^ ' 

The error by which the vector njv(io) violates the periodicity condition ( p3D of the 
n-axis is therefore smaller than 2sec(^/2) tan(^)/(A^+ 1) and converges to for large 
N. 

If one can prove the existence of a suitable number ^ < 7r/2 for some spin transport 
system, one has proven the existence of functions "fiN for this system which satisfy the 
defining equation for the n-axis to arbitrary precision. Since, however, these functions 
fiN do not necessarily converge, this does not prove the existence of an n-axis for such 
a system. 

If the orbit motion can be described in terms of action-angle variables, as is always 
the case for linear motion, and the orbital angle advances for one turn (2-7? times the 
orbit tunes) are not in resonance, then two important conclusions about this tracking 
algorithm with ^ < 7i/2 can be drawn. 

1. If an n-axis n{z) exists, then the sequence njv converges to n{z) linearly in 
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2. If an n-axis exists and the spin rotation angle in one turn is not a linear com- 
bination of orbit phase advances modulo 27r, then the n-axis is unique up to a 
sign. 



The proof is adapted from 0, . The first step will be to show how to define 
a spin rotation angle which depends only on orbital action variables. We assume 
that an n-axis exists and introduce two unit vectors Ui{z) and U2{z) to create a right 
handed coordinate system {n,ui,U2). The vectors ui and U2 are therefore defined up 
to a rotation around the n-axis by an arbitrary phase space dependent angle (f){z). We 
express the spin vectors s in terms of this coordinate system by s = Sifi + S2U1 + 53^2. 
The coefficient Si does not change during the particle motion around the ring since 
the particle transfer matrix R{z) is orthogonal and ensures that {s ■ n) is invariant. 
The spin motion described by the R{z) matrix is therefore simply a rotation around 
the n-axis by a phase space dependent angle i^{z). 




I 1 











cos(z/(2')) sin(z/(z')) 
yO — sin(z/(z)) cos(z/(z')) 



V Si3 



(31) 



If we now introduce the complex quantity s = e '''^*^^)(s2 + ^53) where (\){z) is the 
arbitrary angle, then the spin transport is described by 



■5/2 + iSj-i 



-iu{z) , 



Si2 + isis) 



gi(-i/(i*)+</>(i'))_^. 



(32) 
(33) 



Now we introduce orbital action-angle variables J and $ as well as the angle advances 
Q for one turn around the accelerator. Note that the symbol Q is 2tt times the orbital 
tunes. In these variables the one turn transport is characterized hy Jj = Ji and 
$/ = $j + Q. Since the actions remain invariant during the particle motion, we use 
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— * — * 

the symbols and to indicate the spin rotation angle and the free phase 

of the coordinate system for fixed actions J 

Sf = e^(--/(*)+^/(*)-'^/(«+Q))s. . (34) 

We now show how can be chosen so that the spin motion characterized by 

the exponent becomes simplified. As with any function of phase space, the rotation 
271 periodic in all components $j. Therefore, the rotation angle has 
a periodic contribution and a linear contribution in the phases 

cf>ji^) = ^^j{^)+j.$ (35) 

with some vector j that has integer components. The phase space dependent spin 
rotation e'"^J^'^^ is also a periodic function of phase space. But since on the closed orbit 
(J = 0) the spin motion does not depend on $, therefore z^/($) only has a periodic 
component and no component linear in $. 

If the orbital angle advances Q are not in resonance, then 0j($) can be chosen 
to eliminate the phase dependence of the exponent in equation completely. This 
can be seen by Fourier transformation of the periodic functions z^j($) and 
leading to the requirement 

t>j{k) = IMil - e^'-^) (36) 

for the Fourier coefficients t)j(k) and (p^jik). Therefore, (p^jik) can be chosen to 
eliminate all Fourier coefficients except for A; = 0. With this choice the exponent 
reduces to I'jiO) + j ■ Q- The vector j can in general not be used to simplify this 
expression and we therefore usually choose j = 0. With i/(J) = 2^/(0), the spin 
rotation of equation (0) then simplifies to 

Sf = e-"'^'^s, . (37) 
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We have now achieved the goal of constructing a spin rotation depending only on 
orbital actions. It is interesting to note that if for amplitudes J the integer coefficients 
of j can be chosen so that 



z/(J) + J ■ Q = mod 27r 



(38) 



one can eliminate the spin rotation completely. Here we only analyze the case when 
this resonance condition is not satisfied. 



The stroboscopic average is now performed by the recipe of section [4.1| but this 
time in the coordinate system {n,ui,U2) just constructed. We first estabhsh the 
tracking points cj and note that u = i^icj) is the same for all tracking points. In the 
coordinate system {n,ui,U2), the vector components of the periodic spin no on the 
closed orbit are not constant but depend on the phase space position (J, $). This 
vector is transported from the phase space points Cj to Zq = ($, J) by the rotation 
matrix 

f 1 \ 

(39) 



1 

cos(jz/) sin(jV) 
— sin(jz/) cos(jV) 



leading to the stroboscopic average in this coordinate system 



N 



S/v($, /) 



N + 



1 


i=o \ 





cos(jV) sin(jV) 
- sin(jz/) cos(jV) 



0,1J 



($ - jQ) \ 
(f-jQ) 



0,2.J 



(40) 



The first component of sn is Y.f=o^o iji^ '~ jQ)/i^ + 1)) the second and third 
components in complex notation are 



N 



sn = sn2 + ism 



N + 



-^e-(^%,-.($-jQ), 



(41) 



j=0 



where f{^) = nQ2j(^) +^"'0 3j('^)- terms of the Fourier components f{k) of 



Hq one obtains the inequality 

1 1 _ Q-i{N+l){u+k-Q) 



< 



— rir, r{k)e 

1 ^ 2 



Similarly one obtains for the first component of Sn 

1 1 _ p-i{N+l)k-Q ^ ^ ^ 

MO 



Since |s/v| > cos(^) and ^ < 7r/2, it is guaranteed that sn does not converge to 0. For 
large N, |s/v| is given by ij' "which therefore does not vanish. Since the resonance 
condition (|38|) is not satisfied for any vector of integers j, none of the denominators 
is zero and sn,i and sn converge linearly with 1/N. At this point we impose the 



condition on n that the sum on the right hand side of equation (42) converges. From 
the convergence of Sn the convergence of n^v to n follows from 





1- 1 


n — fiATp = 


\n — 7-^— j- 







\sn 




sn,i 




sn 





< 2 '"",' 7^'^ < 2{p^f . (44) 



\sn 


- Sn,1 




sn,i 





\sn,i\ 

There is a number such that the absolute value of the N dependent part of sn,i 
in equation ( PB| ) is smaller than Uq ^j{0)/2. If the number of turns N is bigger than 
A^*, we conclude 

\nN -n\< 2V2 . '^^j ., (45) 



which, together with equation (^21), establishes that un convergences linearly with 
1/N. 
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The second conclusion to be proven is the uniqueness of the n-axis. In the coor- 
dinate system {n,ui,U2), the periodicity condition ( p3D reads 

1 \ 

cos(z/) sin(z/) ■ n{z) = n{M{z)) (46) 
— sin(z/) cos(z/) / 

with the obvious solution n(£) = (1,0,0)"^ for all z. If another n-axis n2{z) exists, 
then ^2 — n{n ■ 7^2) is nonzero at least at one phase space point and on all iterates of 
this point which can be reached during particle motion. We normalize this difference 
vector at these phase space points and write it as cos(a(z'))'Ui + sin(a(i*))'U2- In orbital 
action-angle variables, the free function dj{^) again has a periodic contribution and 
a linear contribution j ■ $. 

In complex notation, the periodicity condition (^) reads 

Since z^(</) does not depend on $, the periodic part of a has to vanish. The resulting 
a J = j ■ $ can only solve the periodicity condition if there is a vector j with integer 
coefficients fulfilling the resonance equation (|38D and then a second n-axis exists. For 



linear orbit motion with phase advances which are not on resonances of type (|38|) 



equation ( ^Tf ) does not have a solution and the n-axis is unique. There is so far no 



proof for nonlinear orbits, but nevertheless the method of stroboscopic averaging can 
also be used for exploring nonlinear orbit motion. 

4.3 Improved Recipe with Faster Convergence 



In section [4.1| an algorithm was introduced which converged with the square root of 

the number, T, of turns tracked. However, it is possible to obtain convergence linear 

in the T, when one takes advantage of the orthogonality of the spin transfer matrix. 

To illustrate this, we again establish a recipe. 
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1. Define as before 

C = {c,- = {M-^y ■ z,\3 e {0, . . . , TV}} . (48) 

2. Define tliree ortliogonal unit vectors {e^^-*, e*-^-*, e*-^-*}. 

3. Obtain the sets Sj of the three vectors sf \ and by tracking the e^'^) for 
N — j turns 

^ = • . . . • E(cAr) • e-<^) . (49) 

4. Prom the set of vectors Sj and the set Sq, one can obtain the spin transfer matrix 
from the phase space point Cj to Zq denoted by R{cj). This becomes clear when 
one reahzes that — R{cj) ■ sf^ for all k. Obtaining these + 1 transport 
matrices requires the propagation of 3 spins around the circular accelerator for 
N turns. 

5. Now we can compute the set B = {R{cj) ■ no\j e {0, . . . , N}}. This is obviously 
identical to the set denoted by B in the previous section. 

6. The normalized average of B, denoted by fijv, can now again be computed as 
mentioned above and solves the defining equation for an n-axis up to the small 
error discussed previously. 

In this approach one only has to track three initial spin directions over N turns, 
leading to T = 3A^. The error is therefore bounded by 6sec(,^/2) tan(,^)/T. This 
implies convergence linear in the number T of turns tracked. The following example 
illustrates the speed of this method: when the angle ^ happens to be 45°, and we 
require an accuracy at the 10~^ level, this linear convergence approach only requires 
6500 tracking turns; when the angle ^ is small, fewer iterations are needed. 
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4.4 Backward Tracking 



In the two recipes of section ^?T| and mentioned above, we need to find the set C 



of + 1 backwards tracked phase space points, and then launch spins at these points 
and track forward so as to compute the set B. In the case of linear motion, it is trivial 
to obtain these backward tracked phase space points. One simply transforms zq into 
the action-angle variables of the linear motion and determines the phase advance per 
turn of the linear motion. Counting back these phase advances and transforming the 
action-angle variables back into phase space leads to the points Cj. In the nonlinear 
case this would actually require to track for more turns around the ring. 

In the case of the linearly convergent method of section ^.3| , this extra effort 
becomes unnecessary. One can start with the phase space point Zq and launch three 
particles with spins along e'-^^ e^^-*, and e'-'^^ Tracking backward in azimuth defines 
the + 1 sets Pj of the spins Pj^\ ff^\ and ff^^ with 

= R-\cj) ■ . . . ■ R-\ci) ■ e<'^ . (50) 

From the sets Pj and Pq one can again compute the spin transfer matrix R{cj) and 
with these matrices one obtains the set B with the elements bj = R{cj) ■ uq, which 
again leads to ni^ by averaging. 

4.5 Forward Tracking 

There is an even more fundamental problem in the case of nonlinear motion than the 

computation time. When the lattice or the effect of separate nonlinear elements are 

computed by nonlinear transfer maps, the inverses of these maps might not be known. 

In this case the required phase space points Cj cannot be computed at all. Nevertheless 

our method can be used to obtain the vector n^r as follows. The arguments of the 

section on backward tracking can simply be repeated for tracking forward. One 
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— * 

establishes the phase space points Cj = M(cj_i) with cq = Zq for j G {1, . . . , A^} and 
simultaneously the sets Sj by tracking the three unit vectors Sj_i for one turn with 
^) — ^fe) ^Yie fourth step of the recipe in section [4.3| one can then obtain 



the spin transfer matrix R{cj) from the phase space point zq to Cj. The inverse 
of this transfer matrix is simply obtained by transposition leading to the vectors 
= R.{cj)^ ■ no- The normalized average of the vectors bj is then the stroboscopic 
average ninv,N of the inverse motion. For this average the error of the periodicity 
condition for the inverse motion converges linearly to zero. Fortunately, an ra-axis of 
the inverse motion ninv{z) is also an n-axis of the forward motion, since the periodicity 
condition of the inverse motion reads as 

R'\M'\z)) ■ n,nv{z) = ninv{M-\z)) (51) 

and from this it follows that niny{z) also obeys the periodicity condition (^) of the 
forward motion. 

4.6 Tracking Rotation Angles and Rotation Vectors 

One can represent a rotation matrix by its rotation vector (3 and its rotation angle 
G [0, tt]. The simplest representation is in terms of the vector 7 = sin(0/2)/5 and 
K = cos(0/2). The rotation matrix is then given by 

Rij = 2{K,^Sij + 'ji'jj - KSijklk) - Sij . (52) 

This representation can be easily transported through an accelerator by means of the 
equations 

K = K1K2 - 7i ■ 72 , K = \k\ (53) 
7 = (71 '^2 + l2f^i + 72 X 7i)- . (54) 

K 
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It is therefore even more appropriate and faster by about a factor of three to track 
rotation matrices represented by 7 and k directly, than to track three spin vectors. 



5 The Approximation of Linear Spin— Orbit Mo- 
tion 

In this section we consider the special case of linear spin-orbit motion as a way to 



illustrate and confirm the chief features of our algorithm ||T6[. For this purpose we 
define two periodic vectors m and / orthogonal to Uq in order to obtain an orthonormal 
right-handed dreibein at the previously specified azimuth 6q. We then write a spin 
vector s in the form 



'1 - a2 _ + am + pi . (55) 

In the case where both a and (3 are much smaller than unity we may treat the spin 
in a first order approximation in a and (3 so that the spin vector ( ^5)) can be written 

as 

s = no + am + (31 (56) 
and is normalized in linear approximation. We combine the two-component vector 



a 



(57) 



and the orbit vector z into a single vector z of d+2 components. By linearizing the 
equations of motion with respect to z, an initial spin-orbit coordinate Zi at azimuth 
6q is mapped into a final coordinate Zf after one turn around the ring by the 
{d + 2) X [d + 2) spin-orbit transfer matrix M 

Zf = M■z^ , M= ( ^ ^ ) , (58) 



24 



where M again denotes the d x d dimensional orbit transfer matrix. The matrices G 
and Z2 are 2 x d and 2x2 matrices respectively with 

Sf = G-Zi + D-§i. (59) 

In order to simplify the following formulas, we introduce the 2 x d dimensional com- 
ponent of the transfer matrix M for m turns around the ring. We note that 
the two component vector corresponding to uq is no = so that propagating this 
initial spin vector at a phase space point Zi once around the ring leads to the final 
spin Sf = ■ Zi. Therefore, we define the two component vector analogous to 
equation (|26|) by 

N _ / i\/rm ri \ 

M^o) =T.Gm-M-^-z,, M"'=if ^Hn . (60) 

From equation (EHf) it follows that 



m— 1 



G^- M-"" = J2 d!" -G- M~^''+^K (61) 

We expand the phase space coordinate Zq in eigenvectors yj of the one turn transfer 
matrix M 

d 

j=i 

with expansion constants and eigenvalues 

M-f, = A-V, . (63) 

Thus we get 



m—l 

''J 



G^ ■ M-™ ■ = E ^r'^' -G-Vj. (64) 

fc=0 



We do not consider the case of spin-orbit resonances, so that the matrices 1 — 
are nonsingular. Thus equation (|6^ ) simplifies to 



G^ ■ M— ■ = A,(l - A™^-) ■ (1 - \,D)-' -G-y,. (65) 
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Prom this follows 

N 



E G^-M-^-y, = A,[(iV+l)l-(l-Af+iD^+i)-(l-A,D)-i].(l-A,D)-i.G-%- , (66) 



m=0 

SO that 



1 ^ 



Combining this with equations (^Oj) and (|6^) we get 

^(zo) = E - ^J-^)"' -G-Vj , (68) 

where n is the two component vector corresponding to n. For d = Q this yields the 



well known expression of the n-axis for the case of linear spin-orbit motion 0, |T3 
This confirms the method of stroboscopic averaging for the linear case. Also our 
predictions about convergence speed are confirmed as follows. In storage rings the 
particle motion can only be stable if the eigenvalues have |Aj| < 1. To analyze the 
convergence speed, one has to realize that D is a rotation matrix and therefore that 

1(1 - Af+^^^+^)l| < (1 + |A,|^+^)|I| < 2|5| for all 5, (69) 

2 



|n^(io)-n(fo)| < j^Y.\<'Ml-\Dr'-G-y^. (70) 

This inequality shows that also for linear spin-orbit motion one finds convergence 
linear in and it can also be seen that the convergence speed decreases with the 
orbital amplitudes and becomes very slow close to first order spin-orbit resonances. 

6 Numerical Examples 

In order to illustrate which quantities can be computed and how effective stroboscopic 

averaging can be in numerical computations, we apply it to a model accelerator with 

known n-axis. Since this model is somewhat artificial, we also apply our method to 

the vertical motion of the HERA proton ring. For simplicity, the vertical bends at 

the interaction regions have been ignored. 

26 



6.1 Comparison with an Analytically Solvable Model 

In this section we consider a special model with d = 2 and z = {^,J). The 



equations of motion are given by 



-{6) = Q{m,j)xm, (71) 



with 



J) = ^(z/oei + /iv0'(e3sin$ + e2cos$)) , (72) 

where /i is a real parameter. Hence initial coordinates Zi are taken into final coordi- 
nates Zf by Zf = M{zi) with = $i + Q and J/ = Jj. In the following we assume 
1^0 ^ Q and that is not integer. This model corresponds to the rotating field 

approximation often used to discuss spin resonance in solid state physics |T^. We now 
introduce the orthogonal matrix T(e, ip) describing a rotation around a unit vector 
e by an angle (p. Transforming the spin components of s into a rotating frame by 
introducing Srot = H{gi, $) ■ s, one obtains the simplified equation of spin motion 



^{6) = QrotiJ) X Srotie) , QrotiJ) = 7^ \ ^I^/J | • (73) 



If in this frame a spin field is oriented parallel to ^Irot, this field does not change from 
turn to turn. Therefore Urot = ^rot/\^rot\ is an n-axis. In the original frame this 




n-axis IS 



ri($,J) = ^^^_^l((zyo-Q)ei + /iv^(e3sin$ + e2Cos<l>)) , (74) 



A = ^{uo-Qy + fi^J, (75) 

where the 'sign factor' {i^o — Q)/\i'o — Q\ has been chosen so that on the closed orbit 

(J = 0) the n-axis n($,0) coincides with = Ci. As with any function of phase 

space this n-axis is a periodic function in $. 
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We now perform the stroboscopic average by the recipe of section TA_ to compute 
an n-axis at io = ('^', ^) and 60 

1 A' j 



UN 



jJ2I[Ri^~kQ,J)-no (76) 

-L — n 7 — 1 



I N j 
_ 1 

j=0 k=l 

N 



Wrr E n ■ mot, A) ■ T(ei, $) ■ no (77) 



= ^r^2:(ei,-*)-Emot,jA)-ei . (78) 
If A/27r is not an integer one obtains after some tedious manipulations 

|n^(<l>, J) - n($, J)| = V2VT^ , (79) 

^ -n I ^'-^ l-cos((iV+l)A) 

^^-^'+(Ar + l)Vo-Q)^ l-cos(A) ^ • 

One sees that I^at — n| is an oscillating function of J whose local maxima increase 

with J, reflecting the fact that large orbital amplitudes reduce the convergence speed. 

This behavior is plotted in figure |^. In this and the other figures concerned with the 

solvable model we used the parameters uq = O.Gvr, Q = 0.467r, ^ = 0.27r, and $ = 0.32. 

For large equations (|7^) and (^) predict that the convergence is indeed linear in 

as illustrated by the slope of —1 in figure ^. Also one sees that \nN — n\ vanishes 

for A/27r T^integer and if {N + l)A/27r =integer. For A/27r =integer we have 



= (1 + (81) 

Therefore in this case un does not converge to n. This is no surprise since the 
condition A/27r =integer amounts to the resonance condition (^) which leads to this 
non-uniqueness of the n-axis. It is interesting that the local maxima of In^r — n\ do 
not occur at these resonance points. 

Figure ^ shows the variation of the opening angle of the polarization field as 

a function of orbital amplitude. Since uq ■ n = {uq — Q\/ A as well as no ■ un = 
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j^^jLo'^ni'^rou j^) do not depend on the angle variable the depicted angles 
are equivalent to the phase averaged opening angle. The good agreement with the 
analytically computed opening angle of the n-axis shows that accurate values for the 
field can be obtained with a very limited number of turns. 

\n - nioool 



0.0008 




2 4 6 8 10 12 14 



Figure 2: The deviation of the stroboscopic average un from the analytically calcu- 
lated n as a function of the amplitude J in phase space 



This analytically solvable model can also be used to illustrate the construction of a 
phase independent spin rotation angle i^(J). Having got an n-axis, one can transform 
the spin components of s into a rotated coordinate system (n, Ui, U2). With the simple 
choice 

ui{^: J) = ^^(^ X ei) = |''°~^| (e2sin$- egcos^) , (82) 
|nxei| \i^o-Q\ 

■U2($, J) =nxui , (83) 
both of which are clearly periodic in $, one can show that the spin rotation angle is 
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logio(A^) 



Figure 3: Convergence of the stroboscopic average un to the analytically calculated 
n with the number of turns tracked for J = 14. 
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Figure 4: Opening angle of the analytically calculated n (diamonds) and the strobo- 
scopic average fiN for = 20 (+) as a function of phase space amplitude J. 
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independent of $, namely 



However, z/(J) is free up to multiples of the orbit tune as outlined in section |4l^ . We 
can use this freedom to obtain a z^(J) which reduces to uq on the closed orbit (J = 0). 
The only choice which satisfies this condition is obtained by an additional rotation of 
Ui and U2 around n by — 

ui{^, J) = ei^-— sm$ + e2Cos$sm$(--^ — ^ ) 

-^3(1 — -cos <I>H sm^<I>) , (85) 

Wo-Q\ A 

■U2($, J) = fix ui . (86) 

Both are again periodic in $ and again we obtain a spin rotation angle independent 
of viz. 

K^) = r^A + Q- (87) 
I ^0 — V I 

Correspondingly, on the closed orbit ui,U2 reduce to 

«i($,0) = -e3/^^, «2($,0) = 6-2^^5^ . (88) 



6.2 Application to the HERA Proton Ring 

The HERA proton ring stores protons with an energy of 820GeV. At this energy, the 

proton spin rotates about 1566.85 times around the vertical direction in a flat ring 

during one turn. The HERA ring is not completely flat, nevertheless this number 

illustrates the complexity of spin motion. If a change of the phase space position of 

a particle initiates a relative change of the one turn spin rotation by only one part 

in 10^, the angle of precession has changed by 56 degrees. This simple observation 

already hints at the strong dependence of the equilibrium spin direction on the phase 
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space coordinates pointed out in the Introduction. To illustrate this fact, we restricted 
ourselves to vertical motion with a tune of 0.29, ignored the vertical bends to obtain 
a flat ring, and computed the third order expansion of the spin transfer matrix with 
respect to the vertical phase space coordinates using the Differential Algebra code 



COSY INFINITY |T9|. Finally we stored the third order phase space expansion of 
the corresponding rotation vector 7. This power expansion was then used by the 
program SPRINT to obtain the n-axis. Due to the strong phase space dependence of 
spin motion, this expansion does not represent the real spin motion in HERA very 
well, but it is nevertheless useful to demonstrate the applicability of stroboscopic 
averaging. We transported spins by means of this rotation vector and approximated 
the orbit motion by the linear transport matrix. 

In a realistic accelerator, the existence of the n-axis cannot be guaranteed, but 
an approximately invariant spin field can be found, if the series converges. To 
indicate the convergence of this series, we plot |njv — "'^20000 1 for the phase space point 
with Hi = 0.4mm and = in the East interaction region. This corresponds to an 
emittance of GQvrmm mrad and is therefore a particle at approximately Aa of the beam 
distribution. The slope of —1 in the double logarithmic scale of figure |^ illustrates 



clearly that the convergence is linear in as derived in section 

In that section, convergence could only be guaranteed if the angle between the 
vectors 67V and tiq stayed smaller than 7r/2 during tracking. As an example we checked 
this requirement for yi = 0.4mm and found that this condition is violated as illus- 
trated in figure |^, and convergence cannot be guaranteed ad hoc. Nevertheless, an 
approximately periodic spin field is obtained. 

One can easily check this by tracking a spin which is initially parallel to nj^ for 
many turns. In order for n^v to approximate an n-axis, the tracked spins have to lie 
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Figure 5: Indication of Convergence 




Figure 6: Angle between the tracked spins and no for i/i — 0.4mm. To guarantee 
convergence, this angle has to stay between and tt/2. 
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approximately on a closed curve on the unit sphere. Four such closed curves were 
created by computing ni2ooo at the phase space points yi = 0.1mm, yi = 0.2mm, 
yi = 0.3mm and yi = 0.4mm with y^' = and tracking for a further 600 turns. In 
figure we display the projections of these curves on the horizontal plane. 
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Figure 7: The stroboscopic average hn for = 12000 tracked for a further 600 turns 
for yi =0.1mm, 0.2mm, 0.3mm, 0.4mm from top left to right bottom. 

Optimization of the average polarization of a particle beam requires that the 
equilibrium polarization direction for every particle is almost parallel to the average 
polarization direction. We therefore averaged the angle between n^v and no over the 
orbital phases and displayed this divergence for different phase space amplitudes in 
figure |. 

Yet another way of illustrating the importance of the n-axis is illustrated in figure 
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Figure 8: Phase averaged opening angle of the stroboscopic average for = 200 
as a function of the vertical phase space amplitude J. 

^. Particles at 100 different phases at a normalized emittance of about 47rmm mrad, 
corresponding to Ui = 0.1mm and = 0, have been tracked through HERA for 
500 turns while the beam was initially polarized 100% parallel to uq. Similar kinds of 
tracking results have been presented in Since this polarization distribution is not 
the equilibrium distribution, the averaged polarization exhibits a strong beat. This 
figure also shows that when spins at phase space coordinate z are initially parallel 
to n{z), the averaged polarization stays constant. Therefore, by starting simulations 
with spins parallel to the n-axis one can perform a much cleaner analysis of beam 
polarization in accelerators. 

Conclusion 

We have introduced an algorithm for computing the Derbenev-Kondratenko spin 
axis (n-axis) from straight forward spin phase space tracking data with the following 
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Figure 9: Propagation of a beam that is initially completely polarized parallel to no 
leads to a fluctuating average polarization. For another beam that is initially polarized 
parallel to the periodic spin solution n the average polarization stays constant, in this 
case equal to 0.765. 
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features: 



• It can be implemented in any existing spin tracking program. 

• For an accuracy on the 10"'^ level typically less than 3000 turns have to be 
tracked. 

• No artificial damping is needed. 

• Since our method is non-perturbative, no resonance factors appear in the algo- 
rithm and it is applicable even near spin-orbit resonances. 

With the method of stroboscopic averaging important features of accelerators can 
now be analyzed. One very significant field of study will be the computation of 
spin tune spreads, which previously could hardly be analyzed and are now easily 
accessible. Diffusion and damping terms could be added to the analysis and ^dnjd'y 
0] can be computed for the simulation of electron polarization in storage rings. In 
addition numerical checks of uniqueness, convergence near resonances, and the effect 
of nonlinear orbit motion should be made. 

It is the opinion of the authors that now, after the introduction of stroboscopic 
averaging, spin tracking in storage rings should always be initialized with spins parallel 
to the equihbrium distribution, since much clearer analysis becomes possible. 
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